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Synchronization of coupled oscillators is often described using the Kuramoto model. Here we study 
a generalization of the Kuramoto model where oscillators communicate with each other through an 
external medium. This generalized model exhibits interesting new phenomena such as bistability 
between synchornization and incoherence and a qualitatively new form of synchronization where the 
external medium exhibits small-amplitude oscillations. We conclude by discussing the relationship 
of the model to other variations of the Kuramoto model including the Kuramoto model with a 
bimodal frequency distribution and the Millennium Bridge problem. 



I. INTRODUCTION 

Spontaneous synchronization of coupled oscillators oc- 
curs in many biological, physical, and social systems [H- 
0] . The onset of synchronization is often described using 
the Kuramoto model Q where oscillators are directly 
coupled to each other through their phase differences. 
However, in many realistic systems, the coupling be- 
tween oscillators is not direct and instead occurs through 
a shared external medium. Examples include popula- 
tions of synthetically-engineered bacteria Q and single- 
cell eukaryotes communicating through chemical signals 
0, lasers coupled through a central hub 0, and even 
pedestrians walking on a bridge [1, @. In contrast to 
directly coupled oscillators, such systems have received 
relatively little attention [Ifl [ll| . 

In this paper, we present a generalization of 
the Kuaramoto model where oscillators communicate 
through a common external medium. As in the usual 
Kuramoto model, each oscillator is described by a time- 
dependent phase, 9i[t), which in the absence of coupling, 
rotates at its natural frequency LOi. We concentrate on 
the case when the number of oscillators is large, ^ 1 
and the natural frequencies are assumed to be drawn 
from a prescribed distribution function g{uj) with mean 
frequency ojq. The oscillators are coupled through an ex- 
ternal medium which has an amplitude and phase and 
is described using a complex number Z{t) = i?(i)e'*'*'. 
The model can be derived from the system introduced in 
[Tl| to study dynamical quorum sensing by considering 
the weak interaction limit where the amplitude of each 
individual oscillator's limit cycle remains approximately 
constant. 

The dynamics of the external medium gives rise to in- 
teresting new phenomena not seen in directly coupled 
oscillators. A key dynamical parameter is the density, p. 



of oscillators. For large densities, p — >■ cxj, our model 
reduces to the usual Kuramoto model. Surprisingly, 
in the limit p — > 0, the phase diagram of our model 
can be mapped to the phase diagram of a Kuramoto 
model with a bimodal frequency distribution (l^ with 
the added restriction that the dynamics is constrained 
to lie within the Ott-Antonsen manifold [l^. At low 
but non-zero densities, p ^ 1, the system exhibits bista- 
bility between incoherence and synchronization as well 
as between two difference kinds of synchronized states. 
Additionally, when wq = 0, the dynamics of our model 
is in one-to-one correspondence with that of the Millen- 
nium Bridge problem. Thus, the Kuramoto model with 
coupling through an external medium represents a par- 
ticularly simple physical model which captures behavior 
exhibited by a large number of other Kuramoto-like mod- 
els . 

The outline of the paper is as follows. In section |lll 
we define our model and show how it naturally arises 
from considering the weak-coupling limit of the dynamic 
quorum sensing model considered in In section Hill 
we use the Ott-Antonsen ansatz to derive equations for 
the dynamics of our model for an arbitrary oscillator- 
frequency distribution. In the next section IIVI we con- 
sider the special case of Lorentzian distribution where 
we can derive analytic equations for the stability of var- 
ious phases. In section |Vl we use these equations and 
numerical simulations to construct a phase diagram for 
our model. In section IVI[ we discuss the relationship 
of our model to other problems including the Kuramoto 
model with bimodal frequency distribution and the Mil- 
lennium Bridge problem. We then discuss the implica- 
tions of model and conclude. 



II. DERIVATION OF THE MODEL 
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The Kuramoto model with oscillators coupled through 
a shared external medium can be derived by considering 
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the weak coupling limit of the model introduced in 
The model in (DQS) consists of a population of limit- 
cycle oscillators, ^i, each with a natural frequency ujj, 
diffusively coupled to a shared medium, Z . When chem- 
icals leave the oscillators and enter the medium, they are 
diluted by a factor a = Vint/Vext <C 1, the ratio of the 
volume of the entire system to the that of an individual 
oscillator. Furthermore, the external media is degraded 
at a rate J. The dynamics of the system arc captured by 
the equation 



dzj , , 

— = (Ao + luj, 



Diz, 



Z) 



Z) - JZ 



(1) 

(2) 



z is one when all the oscillators have the same phase and 
zero when they are completely out of phase. Notice that 
the external medium communicates with the oscillators 
only through z = re"^ and we can rewrite ([7|) in terms of 
the order parameter as 



dZ 



= pD{z- Z) - {J + iojQ)Z 



(9) 



III. DYNAMICAL EQUATIONS FOR 
ARBITRARY FREQUENCY DISTRIBUTIONS 

A. Thermodynamic limit and the Ott-Antonsen 
Ansatz 



with Ljj drawn from a distribution h{uj) which we assume 
to be an even function about a mean frequency wq- By 
introducing a dimensionless density, p = aN, and shift- 
ing to a frame rotating with frequency ujq, we can rewrite 
(21) above as 



dZ pu , 



pD 

N 



Z) - {J + iujo)Z, 



(3) 



where the frequencies Wj are now drawn from a distribu- 
tion g(uj) with mean zero. 

We are interested in the limit Aq ^ Z? where individual 
oscillators are weakly coupled. In this limit, , the ampli- 
tude of the limit cycle is not modified by the coupling and 
the dynamics is well-described by modeling each oscilla- 
tor by a single phase variable [IJ]. Explicitly, rewriting 
^ in polar coordinates with Zi = Vie^^' and Z = i?e** 
yields 



±1 
dt 
dO^ 

dt 



{Xo-D-r^)rj+DRcos{<S>-ej) (4) 
DR 

H sin($-6'j). (5) 



In the large Aq limit, the first term dominates the right- 

1 /2 

hand side of ^ and thus ~ Aq in steady-state. Dcfin- 



,,1/2 

mg n = r-i/Ao' 



1 /2 

1 and R = R/Xq , we are left with 
the reduced equations that define our model: 



d0i 

-^=u,+DRsm{<^-e,) 



dt 



N 



(6) 
(7) 



We will drop the tilde in the remainder of this paper and 
note that r = 1 will correspond to the fully synchronized 
state. Interestingly, however, we will find that 

Finally, it is useful to define the Kuramoto order pa- 
rameter, 



N Z^ 



(8) 



In this and the next section, we derive equations for 
the time-dependence of the distribution function as well 
as the steady-state values of the order parameter within 
the low-dimensional manifold of states introduced by Ott 
and Antonsen [l^. We have explicitly checked that the 
dynamics of our model is captured by the Ott-Antonsen 
ansatz using numerical simulations. In what follows, we 
restrict ourselves to the thermodynamic limit, N oo. 
In this case, the dynamics is well described by the time- 
dependent density function f{9,Lu,t), with the fraction 
of oscillators of frequency w lying between 9 and 6 + d9 
given by fdO. 

The evolution for / in time is governed by the conti- 
nuity equation 



(10) 



where dot signifies a derivative with respect to time. 
Plugging in (|6]) and rewriting the result in complex no- 
tation gives 



Z*e''')f 







(11) 



Note that unlike the ordinary Kuramoto model the 
external medium, Z, enters in the second term in place 
of the continuum limit of the usual order parameter, z = 
J /e*^. Within the Ott-Antonsen ansatz, the dynamics 
lies on a submanifold where the distribution functions are 
of the form, 



f{d,io,t) 



.g(^) 
27r 



(12) 



where a{uj,t) is a function of lu and t, g{Lo) is the fre- 
quency distribution of the oscillator population, and c.c. 
denotes complex conjugate. Plugging this ansatz into the 
continuity equation gives 



da 
'dt 



D 



■ loja 



{Za- - Z*) 







(13) 



The Ott-Antonsen ansatz can also be used to derive 
a simple expression for the order parameter z. In the 
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thermodynamic limit, we can write the order parameter 
as 



z{t) = / dLjdefie,uj,t)e" 



(14) 



Plugging in ()12|) and noting that only the term propor- 
tional to e^*^ is non-zero yields the relation 



z = y dujg{uj)a* . (15) 
We will use this definition extensively in what follows. 

B. Stability of the incoherent state 

We start by examining the stability of the incoherent 
state. The calculation is a straightforward generaliza- 
tion of that performed on the original Kuramoto model 
p^ . In the incoherent state, the phase of the oscilla- 
tors is uniform and / = Notice that corresponds to 
choosing a{u;,t) = in ijl2|) . This choice satisfies the 
dynamical equations ([T3|) . ([9|), and ([T5| when Z = z = 0. 
To calculate the stability of this state, we linearize the 
dynamical equations about the incoherent state by sub- 
stituting / = 1/2tt + efi into ((T0| and keep terms linear 
in e. This is equivalent to keeping only terms linear in a 

in (nr 



da D 



0. 



(16) 



For stability analysis, it is sufficient to consider expo- 
nential solutions of the form a = Acxp[{X + in)t] and 
Z = Zoexp[(A — with A and real numbers. The 

parameter A determines the stability of the incoherent 
state and is a rotation frequency with > 0. In par- 
ticular, in the lab frame, the steady-state solutions rotate 
with a frequency wq — slower than the mean frequency 
Wo due to an external medium induced "drag" . This is 
in contrast with the usual Kuramoto model with direct 
coupling where £7 = 0. Plugging in solutions of this form 
into yields 



A 



DZ* 



2(A + + iuj). 



(17) 



Using this in eqns. pSj) and © results in a complex 
equation for Q, and A of the form 



A + i(a;o - ^) + pD + J 



dujg{uj) 



(18) 



To determine the boundary of stability, we put A = 0+ 
in the equation above and equate the real and imaginary 
parts to get 



pD + J 

+ Wo 

pD^ 



T:g{n) 



P 



du 



.g(^) 

u! + n 



(19) 
(20) 



In arriving at this equation, we have used the standard 
identity = P{-^) — iTrS{x) when e = 0+, with P 

denoting principal value. Equation ([20)) gives implicit 
equations for the surface in the p-D-J parameter space 
where the incoherence state changes stability. We will 
use this result to construct the phase diagram for our 
model below. 



C. General Locked Solution 

For general g(w), we derive an equation governing 
locked states as follows. For such solutions, we expect 
that the order parameter and external medium rotate 
at a constant frequency — H, or equivalently loq — 
in the lab frame. Wc look for solutions of the form 
a = A{uj)e^^^ , Z ~ Zoe~*^*, and z = zoe~'^^*. Plugging 
these functional forms into Eq. (fT51) gives 



Defining b ^ oj + fl and R ~ \Zo\ we see that 



Aib) 



ib±y/-b^ + D^R^ 



(21) 



(22) 



To ensure the dynamics stays on the Ott-Antonsen 
manifold, its necessary to make sure that |A(a;)| < 1 
for all io. In particular, this inequality ensures that the 
sum in ([T2|) does not diverge. This can be guaranteed 
by properly choosing the positive or negative solutions 
to (|22p as follows. Defining 



A = x/DZo 



it suffices to choose 



X{b) 



ib + V 
ib±V 



'-b^- 


^D'^R^ 


'-b^- 


^D'^R^ 



D^R^ 



(23) 



b> DR 

-DR <b<DR (24) 
b < -DR 



-ib-V^ 

Finally, we can plug in this ansatz into ^ to get 

[i{ujo-n) + pD + J]R^ = DZop J db'g{b'-n)A{b') (25) 

Together, ([M]) . and (^5)) give a dispersion relation- 

ship between and R for locked solutions. 



IV. THE LORENTZIAN DISTRIBUTION 

A. Dynamical equations 

In general, the general equations (|25p derived in the 
last section must be solved numerically to determine the 
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FIG. 1. Phase diagram as a function of parameters. A = 0.1 
in all figmes. (top) ujo = 0.0577. (middle) ujo = 0.4 (bot- 
tom) uq — 10. M stands for monostability and B stands for 
bistability. Ml - incoherence, M2 - coherent oscillations, B - 
bistability between incoherence and coherent oscillations, and 
B - bistability between large amplitude phase-locked solution 
and small amplitude phase-locked solution. Note the change 
of scale and location of both the horizontal and vertical axes. 



existence of locked solution. However, just as in the ordi- 



nary Kuramoto model, the mathematic simplifies when 
f;(cj) is a Lorentzian distribution 



gi^) = - 



1 



TT ^2 -f A2 ■ 



(26) 



For Lorentzians, the integrals over giiS) appearing in 
eqns. (|T9)) . (|20l) . and (|25|) can be performed explicitly 
and one is left with algebraic equations These algebraic 
equations can then be solved to yield the phase bound- 
ary of the incoherent state and the existence of locked 
solutions parameterized by the frequency and amplitude 
of the external medium, {fl,R). One subtlety that re- 
mains is that we are not able to analytically determine 
the stability of locked solutions (except for the incoherent 
solution), and we therefore resort to numerical solutions 
to test stability. 

The incoherence boundary can be computed explicitly 
by plugging in ([^ into ([T^ and using contour integra- 
tion. This yields. 



A + pD + J 



A + pD + J 



1 - -. \ +A{pD+J)~pD^/2 ^ 

(27) 

The dynamical equations for locked solutions can also be 
calculated explicitly for a Lorentzian distribution. The 
key simplification is that the integral in (|15p can again 
be performed by contour integration so that 



z(t) = a*{-iA,t). 



(28) 



Thus, evaluating eq. ([13)) at w = — lA, we obtain the fol- 
lowing equation for the dynamics of the order parameter 



dz . D , ^ ^ 



Z) = 0.. 



(29) 



When combined with the equation for the external 
medium (|9]), these equations completely specify the dy- 
namics of the system on the Ott-Antonsen manifold. 

To look for locked solutions, we plug in the functional 
form z = re'^* into (|9]) . This implies that 



Z{t) 



pDr 



i{n + ujo) + pD + J 



(30) 



Inserting the solution for Z{t) back into Eq. gives a 
complex equation for r and the rotation frequency $7, 



J 



{in + A)r- 



rpD 



2 r 



pD + J - i{VL + uja) pD + J + i{n + uJo) 



(31) 



Locked solutions are solutions to this equation with r > 
0. Notice, that even though the incoherent state, r = 0, 
is always a solution, it can be stable or unstable, and the 
stability boundary is given by the phase boundary (j27p . 



We now look for other locked solutions, i.e. those with 
r > 0. To do so, we divide ([3T|) out by r, and solve for 



5 



as a function of g = to get 



n{q) = -uJo + pD\l^{l-q) 



1. 



(32) 



Substituting the solution back into pil) yields a cubic 
equation for q = . In general, this expression becomes 
unwieldy and must be solved numerically. However, we 
can still draw many conclusions. Physical solutions must 
have real values of 51 and real, positive values ior q > 0. 
Since the end result is a cubic equation for q, there can be 
0,1,2, or 3 locked solutions for a given parameter range. 
Each of these solutions can be stable or unstable. We 
will use these facts to construct a phase diagram in the 
next section. 



B. Equations in terms of phase difference 

It is useful to rewrite the dynamics of our model in 
terms of the phase difference, ip, between the order pa- 
rameter and external medium. Let Z = i?e'*, z = re^'^, 
and ip = ^ — (j). In terms of these variables, we can rewrite 
the dynamical equations of motion ((^ and ©, as, 

D 



r = -Ar + —(1 - r^)Rcos'ip 

R = -{pD + J) R + pDr cos 
R 

■0 = -Wo - D— sin-iA (2p + + 1) 
2r ^ 



(33) 



These equations fully characterize the dynamics within 
the Ott-Antonsen manifold. 

For the special case of locked solutions, there is a rela- 
tionship between the phase difference a-nd the rotation 
frequency J7. In particular, by examining the phases of 
we have 



cos-0 



pD + J 



yj[pD + JY + {^i + uoy- 



(34) 



Alternatively, taking the absolute value of ((30|) yields a 
relationship between the magnitude of the order param- 
eter, r, and the amplitude of the oscillations in the ex- 
ternal medium, i?. 



R 



pD 



pD 



^{pD + JY + {n + uj(>Y pD + j 



rcos^. (35) 



Notice that this is consistent with since R ~ for 
locked solutions. 



V. PHASE DIAGRAM AND BISTABILITY 

A. Constructing phase diagram from dynamical 
equations 

In this section, we discuss the phase diagram for our 
model. The dynamics of the external medium leads to a 



much richer phase diagram than that for the Kuramoto 
model with direct coupling. For brevity, we focus our 
discussion on the case where g{uj) is a Lorentzian distri- 
bution. We have numerically verified that a similar phase 
diagram occurs for other distributions. Furthermore, to 
reduce the number of parameter, we concentrate on the 
case where A = 0.1 and J = 0. Note that, in contrast 
to directly coupled oscillators, the solutions exhibit non- 
trivial A dependence because A must be compared to 
the mean frequency ujq- On the other hand, the behav- 
ior for J 7^ is expected to be qualitatively similar to 
the J = case because the steady-state solutions in eq. 
([3T|) depend on p, D, and J only through the two control 
parameters pD + J and pD^ . 

We can use the results from the previous sections to 
derive the phase diagram in the (p-D) plane for various 
values of loq (see Fig. [1]). The red line is the incoher- 
ence stability boundary given by ^7} . The incoherent 
solution is stable to the left of this line and unstable to 
the right. Recall, that in addition to the incoherent solu- 
tion, there can exist 0,1,2, or 3 locked solutions which can 
be stable or unstable. These locked solutions are given 
by the roots of cubic equation for q ^ resulting from 
substituting ()32p into (|3T|) . Furthermore, notice that the 
only way that the number of locked solutions can change 
is when one of the roots of the corresponding cubic equa- 
tion changes sign or if new real ones appear. In order for 
the former to occur, the solution corresponding to the 
root changing sign should collide with the incoherent so- 
lution r = and change its stability. Thus, this can occur 
only at the incoherence boundary given by the red line. 
This implies transitions at the incoherence boundary are 
continuous. On the other hand, new real solutions appear 
when the discriminant of the corresponding cubic equa- 
tion equals zero. Note that the phase boundary defined 
by discriminant equalling is the only place where discon- 
tinuous phase transitions can occur. That is, real positive 
solutions can come into existence here at non-zero am- 
plitude. This occurs, for instance, in the transition from 
Region Ml to B, in the middle panel of Fig. ([1]). These 
lines are plotted in blue in the phase diagrams. The num- 
ber and type of stable equilibria must be identical within 
each of the regions carved out by the red and blue lines, 
allowing for a straight forward construction of the phase 
diagram. 

To check these arguments, we ran numerical simula- 
tions of = 1000 oscillators with frequencies drawn from 
a Lorentzian distribution and measured the order param- 
eter after decay of initial transients. We initialized the 
system with random phases and chose the initial magni- 
tude of R to be cither zero or a finite value to probe the 
stability of the incoherent or partially locked state, re- 
spectively. The results were in excellent agreement with 
the phase diagram in Fig. [1] We also studied a Gaus- 
sian .g(cj) distribution and similarly found that the locked 
solutions exist within the low-dimensional manifold and 
the system exhibits qualitatively similar behavior to that 
found for a Lorentzian distribution. 



B. Phase diagram and bistability 

We now discuss the phase diagram in greater detaih 
The top panel in Figure [T] shows that the typical phase 
diagram for small loq- The phase diagram is similar to 
that for a directly-coupled Kuramoto model. There are 
two phases, an incoherence phase and locked phase, with 
a density-dependent critical coupling D that marks the 
transition between the two phases. As usual, the order 
parameter r is zero for the incoherent phase and close to 
one for the locked solution. 

When cjQ is increased, shown in the middle panel of 
Fig. [1] a bistable region between incoherence and coher- 
ence appears at low densities. This region results from 
the subtle interplay between the "inertia" of the exter- 
nal medium and the amplitude of the order parameter. 
The influence of the oscillators on the external medium 
is proportional to the density and the amplitude r of the 
order parameter; see ([M]). Thus, if the oscillators are 
incoherent they cannot entrain the media. For large r 
the opposite is true, giving rise to the bistable region. 
At higher densities the phase diagram is topologically 
similar to the ordinary Kuramoto model with a direct 
transition between incoherence and coherence. 

The phase diagram develops more features as ujq in- 
creased further (bottom of Fig. [ij . In addition to all of 
the behaviors discussed above, there also exists a region 
of bistabiity between two different locked solutions: one 
where the amplitude of the external medium oscillations 
is small and another where the amplitude is large. In 
both cases, the amplitude of the order parameter r is 
close to one. The amplitude R and the phase diflFerence 
between the order parameter and external medium, t/i, 
for these locked solutions can be calculated directly us- 
ing the aforemetioned cubic equation and are shown in 
Fig. [21 Notice that the external medium oscillates out 
of phase with the order parameter for the low amplitude 
locked solution. Finally, we note that the size of this 
bistable region increase as coq is increased further. 
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FIG. 2. Amplitude, R, of the external medium (top) and 
magnitude of the order parameter, r, (bottom) for non-trivial 
locked solutions as a function of coupling D, for ojq = 10, 
and A — 0.1. The black lines represent stable solutions and 
gray lines unstable solutions. Notice the bistable region cor- 
responding to region B in Fig. [T] Insets show the locked os- 
cillation frequency in the rotating frame (top) and the phase 
difference, ^p, between the medium and order parameter, z. 



can also be seen in the inset of Fig. [2] 



C. The low amplitude locked solution 

A novel aspect of coupling oscillators through an exter- 
nal medium is the existence of the aforementioned locked 
solution where the order parameter amplitude r is nearly 
unity but the amplitude of the oscillations in the exter- 
nal medium R is small. Interestingly, this locked solu- 
tion always appears together with the usual locked solu- 
tions where both r and R arc close to one. Figure ?? 
shows numerical simulations confirming the analytic pre- 
dictions from the Ott-Antonsen ansatz. This new type of 
locked solution is possible because the effective coupling 
between oscillators becomes small in ^ when R is small 
even when D is large, allowing oscillators to rotate near 
their natural frequency. Furthermore, pSI) implies that 
the order parameter and external medium will oscillate 
at nearly 7t/2 out of phase for such a locked state. This 



VI. RELATION TO OTHER SYSTEMS 

The model studied here is closely related to other vari- 
ants of the Kuramoto model. For ujq = 0, the steady- 
state dynamics of our mode is identical to that of the 
Millenium Brigdge Problem 0, Q . Somewhat more sur- 
prisingly, in the p — > limit, the steady-state dynamics 
of our model can be mapped onto the dynamics of the 
Kuramoto model with a bimodal frequency distribution 
for the special case when the dynamics of the latter 
model are restricted to lie on the Ott-Antonsen manifold. 
We discuss both of these mappings below. 

A. Millenium Bridge Problem 

To understand wobbly behavior of the millenium 
bridge, the authors of represented the dynamics of 
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pedestrians walking on a bridge using a simple mathe- 
matical model in which the bridge is modeled as a driven 
harmonic oscillator. Pedestrians are also modeled as os- 
cillators that try to phase lock with the bridge. We now 
show that the steady-state, locked solutions of the class 
of systems to which the Millenium Bridge belongs can be 
mapped onto the solutions of our model. 

In the Millenium Bridge Problem, the dynamics of the 
bridge position, X, is governed by the equation 

M— + B—+KX = GY^sme, (36) 

i=l 

with M, B, and K the mass, damping, and stiffness of 
the bridge, respectively. The pedestrian's footsteps, rep- 
resented by the 9i. It is also useful to write X = i?sin$. 
Notice, that for locked solutions where R = est. and $ = 
fit, we can relate the bridge position in the Millenium 
Bridge Problem to the external medium in the Kuramoto 
model by noting that formally we can write X = lm{Z}. 
Furthermore, notice that in the non-rotating, lab frame, 
isolating the imaginary part of Eq. ^ gives the equation 
m cos $ + (pL* + J)i? sin $ = ^ sin Oj . 

Comparison with the equation resulting from plug- 
ging in rotating solution X = i?sin$ into (p6| yields 
a mapping between the steady-state dynamics of the two 
problems. In particular, it is clear there exist a map- 
ping between parameters (p, £>, J, wq) — >■ {G,B,K,M) 
for which R and il are preserved. In particular, assum- 
ing i? 7^ gives the mapping pD = aG, B = 1/a, and 
pD + J = a{K — Mft^), where il{p,D, J,u!o) is a function 
of parameters and a ^ is an arbitrary constant. It is 
interesting that, due to the freedom in specifying a, the 
bridge can be under-, critically-, or over-damped. 

In the original Millenium Bridge Problem, the authors 
restricted their considerations to case where the natu- 
ral frequency of the pedestrian and bridge are identical 
0. The rational for this was the "worst case scenario" 
for wobbling of the bridge. Mathematically, this corre- 
sponds to choosing = in the equations above. For 
this choice, the phase diagram is identical to that of the 
ordinary Kuramoto model with incoherent and coherent 
phases (see top of Fig. [IJ and this was what was found 
in Q. However, as discussed in the last section, the more 
realistic case where the natural frequency of the pedes- 
trians and bridge differ so that 7^ 0, gives rise to a 
much richer phase diagram seen in this problem. 

B. Bimodal distribution 

Another variant of the Kuramoto model studied re- 
cently using the Ott-Antonsen ansatz is the Kuramoto 
model with a bimodal frequency distribution. The Ott- 
Antonsen dynamics captured the transition from coher- 
ent to incoherent states, including bistability for certain 
parameters, but failed capture the full dynamics richness 
of the bimodal model. For example, the Ott-Antonsen 



dynamics did not allow for standing wave solutions which 
were seen numericall y. W ithin the Ott-Antonsen mani- 
fold, it was argued in [121 that the dynamics of the system 
is well described by the set of equations for the square 
of the magnitude of the order parameter g = and a 
phase difference tp between the phases of the oscillators 
locked around zbwo (see [T^ l 

q= -Aq+{q-q'^)[coai} + l]] 
= cJo - (1 + g)sin-!/i, (37) 

where the dot indicates a derivative with respect to t = 
D/2t, A = (4A)/L>, and wq = (4a;o)/£'. 

We now show that the steady-states of the model con- 
sidered in this paper are identical to the those derived 
from the equations above. To do so, it is useful to work 
with the dynamical equations ([55)1 in terms of R = \Z\, 
r = |z|, and the phase difference -0. Furthermore, we con- 
centrate on the case when J = 0. Since we are interested 
in steady-states, we can set 7? = in ([33]) and derive 
a relationship between R and r given by i? = rcos'0- 
Substituting this into the remaining two equations gives 

= 7' = -Ar + ^ r^)r cos^ 

= i' = -ujo - y tan t/i (2p + r^i + 1) . (38) 

Writing q = r'^ and taking the limit p — > gives 
= q = -2Aq + Dq(l - q) cos^ ip 
= = -Wo - jsin2V'(g + 1). (39) 

Finally, making the substitution ijj = 2-0, t ~ D/2t, A = 
{AA)/D, and uJq = (4ajo)/-D yields the equations ([57)1 . 

The mapping can also be derived by directly examining 
equations ([27]) and ([5T|) . In the limit p — > 0^+, the 
stability boundary of the incoherent state ((27)) reduces 
to 

i^ = 2(A2+cj2)/A. (40) 

This is identical to the boundary calculated for directly- 
coupled oscillators with a bimodal frequency distribution 
in [12| . We can also directly derive an equation for the 
amplitude of frequency-locked solutions for p = 0^ by 
substituting the (|32)) into ([3T|)add expanding to lowest 
order in p. This procedure yields the equation 



u;o = ±^^^^iD-Dq-2A). (41) 

As expected, this is identical to the amplitude equa- 
tion for directly coupled oscillators with a bimodal fre- 
quency distribution under the identification w = 4a;/D, 
A = AA/D. Thus just as in the bimodal case, multiple 
solutions can arise in our model via a saddle-node bifur- 
cation which occurs when duJo/dq = 0, or equivalently 
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FIG. 3. Relationship of the Kuramoto model with external medium to other variations of the Kuramoto model. 



q ~ 2 — ^ — -. Plugging this expression into (PT|) 

gives the location of the saddle-node bifurcation with the 
caveat that when the bifurcation occurs to the left of the 
incoherence stability boundary, the resulting q is nega- 
tive and therefore unphysical. By equating this curve 
with Eq. (|41]) . we find new physical solutions with posi- 
tive q arise when wo > A/ -s/S. 

Somewhat surprisingly, we have shown that the steady- 
states of our model in the p, J — s> limit is identical, 
with the same parameters, to that of a Kuramoto model 
with bimodal distribution restricted to the Ott-Antonsen 
manifold under the identification of the phase difference 
between the order parameter and external medium in our 
model with the twice the phase difference between the 
order parameters of oscillators locked around iwo in the 
bimodal model. This hints that the essential physics of 
both problems is contained in two coupled oscillators. It 
would be interesting to explore this further in the future. 

VII. DISCUSSION 

In this paper, we have studied the a variation of 
the Kuramoto model where phase-oscillators are coupled 
through a common external medium. Such a model is 
likely to be widely applicable to biological and phys- 
ical systems where oscillators communicate with each 
through some chemical signal or physical membrane. An 
important distinction between the model studied here 
and the Kuramoto model is that there is an important 
new control parameter, the density of oscillators in the 
medium. This allows for interesting new effects such as 
a density-dependent transition to oscillations which has 
been dubbed Dynamic Quorum Sensing [T6l.[T7j. We used 
the Ott-Antonsen ansatz in combination with numerical 
simulations to investigate the dynamics of this model. 



We found that the model had a rich phase diagram with 
bistability between incoherence and locked solutions, as 
well as between two different types of locked solutions. In 
addition, as summarized in Fig. [31 the model is closely 
related to other variants of the Kuramoto model in vari- 
ous limits. 

The underlying reason for the complex phase diagram 
in our model is the dynamics of the external medium. 
The external medium has two distinct effects. First, it 
introduces an effective time-delay for communication be- 
tween oscillators. This time delay manifests itself mathe- 
matically by noting that 51 7^ 0. It was previously shown 
that introducing a fixed time delay leads to bistability be- 
tween different locked solutions as is found in this model 
[isj . However, since the time-delay is not fixed in our 
model, we do not see the hierarchy of locked solutions 
seen in the direct coupling model with delay. Second, 
the medium has an "inertia" so that the a natural fre- 
quency becomes important. In particular, it is precisely 
when the natural frequency of the external medium is 
large that the phase diagram of our model differs the 
most from the ordinary Kuramoto model. 

The dynamics of the model presented here are equiv- 
alent to the Millenium Bridge problem when ujq = 0. 
This is unsurprising since it the bridge acts as an effec- 
tive medium through which walkers communicate. What 
is somewhat unexpected is the richness in the dynamics 
that emerges when the resonant frequency of the walkers 
and bridge differ (i.e. ujq 7^ 0). More surprisingly, the 
model with an external medium interpolates between a 
directly coupled Kuramoto model with unimodal and bi- 
modal distributions as a function of the density p, with 
p — >■ corresponding to a bimodal distribution, and 
p — > 00 a unimodal distribution (see Fig. An impor- 
tant caveat is that this is only true when the dynamics 
of the Kuramoto model with bimodal distribution is re- 
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stricted to lie on the Ott-Antonsen manifold. This hints 
that the underlying the dynamics of many Kuramoto 
models is captured by the Ott-Antonsen ansatz because 
the steady-state dynamics reduces to that of a few cou- 
pled oscillators. The failure of the ansatz to capture the 
standing waves in the Kuramoto model with bimodal dis- 
tribution is likely because this violates the simple picture 
above. It would be interesting to understand the connec- 
tion between various Kuramoto-like models further. 

Perhaps the most experimentally interesting finding of 
the paper is the predicted bistability at low densities. It 
would be interesting to see if this could be observed in 
an experimental system. The largest obstacle to this is 
that in systems where oscillators have both a phase and 
amplitude, bistability is masked by an oscillator death 
phase where the amplitude of all oscillators is pulled to 
zero pj| . Thus, any experimental realization would re- 
quire that the amplitude of oscillators remain be held 
fixed and the phase oscillator approximation apply even 



at strong couplings. Thus, it is unlikely that bistabil- 
ity exists in experimental setups similar to those used to 
study dynamical quorum sensing such as the BZ reac- 
tion with catalytic beads p^j and quorum-sensing cou- 
pled bacteria [5[. Nonetheless, it will interesting to see if 
the results here can be experimentally tested. 
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